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^ . Abstract 
(N 

The purpose of this paper is twofold. An immediate practical use of the presented 
algorithm is its applicability to the parametric solution of underdetermined linear 
ordinary differential equations (ODEs) with coefficients that are arbitrary analytic 
functions in the independent variable. A second conceptual aim is to present an 
algorithm that is in some sense dual to the fundamental Euclids algorithm, and thus 
an alternative to the special case of a Grobner basis algorithm as it is used for solving 
linear ODE-systems. In the paper Euclids algorithm and the new 'dual version' are 
[~». ' compared and their complementary strengths are analysed on the task of solving 

underdetermined ODEs. An implementation of the described algorithm is interactively 
accessible under [7]. 
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1 Introduction 



Underdetermined ordinary and partial differential equations (ODEs/PDEs) are typical ob- 
jects of investigation in control theory but also overdetermined systems of equations where 
originally the number of equations is larger than the number of functions may turn into 
underdetermined problems in the course of their solution. Examples are the conditions for 
Lie-symmetries and conservation laws of linearizable differential equations. These overdeter- 
mined systems of conditions have solutions that involve arbitrary functions of one or more 
variables, i.e. in the process of solving these systems towards its conclusion, underdetermined 
problems often occur. 

Algorithms for solving underdetermined linear equations and systems (ODEs, PDEs, 
multidimensional discrete,..) are known ([1],[2],[3]). There also exists efficient software 
by Daniel Robertz et al. [I], [5]. The purpose of this paper is to describe an alternative 
algorithm that is elegant, efficient and in some sense complementary to the fundamental 
Euclids algorithm which is the basis for the Grobner basis algorithms that are used in other 
implementations. An open question is whether the new algorithm can be generalized to 
partial differential equations (PDEs) and thus not only be a complimentary algorithm for 
the non-commutative differential Grobner basis algorithm for solving ODEs but also for 
solving PDEs. 

An implementation of algorithms described in this paper is accessible online under [7] 
and is part of the package Crack ([5]) for solving overdetermined systems. As we will 
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describe, a new feature of these programs are flags which allow to prevent denominators or 
that allow to reduce the number of terms in the solution. This can be achieved by absorbing 
explicit x-dependent factors into new functions that are introduced during the computation. 
We consider a single linear ODE 

= f^DJi + c 00 (x), r>l, (1) 

i 

A = L w) I ^ J ' Ui > 



3=0 



for functions fx ■ ■ ■ f r of the independent variable x. The equation may be homogeneous 
(coo = 0) or inhomogeneous (coo ^ 0). With the requirement that the ODE is underdeter- 
mined we only assume that at least two functions fi are involved, i.e. r > 1. The differential 
order of each f is nj. For simplicity we assume that all coefficients are sufficiently often 
differentiate, at least X)i n % times. An upper index in round brackets indicates the number 
of differentiations, i.e. ^£ = fi and low order derivatives are denoted by apostrophe, e.g. 

d 2 fi _ ffl 
dx 2 Ji ' 

The task is to find the general solution of the ODE in the form of explicit differential 
expressions Fj 

fi = F i (x,h 1 {x),...,h r ^ 1 {x)), z = l,...,r (2) 

in terms of parametric functions h\{x), . . . , h r _i(x) which are either all free or one of them, 
say hi(x), having to satisfy an ODE and the other hj(x) being free. 
For example, the general solution of the ODE 

= fx 2 + g"x - g'x 2 + f + 3x 

for / = f(x), g = g(x) can be written in the form 

/ = -= = — A = ((x 5 -x 3 + 3x)h" -(x 6 + x 4 + 3x 2 -Q)ti (3) 

J x s - 2x e + 7x 4 - Qx 2 + 9 V v ; v 

+ (3x 5 + 3x 3 + I7x)h - 3x 8 + 3x 6 - 16x 4 + 9x 2> 
x 

9 ~ 2(x 8 - 2a; 6 + 7a; 4 - 6a; 2 + 9) 

+ (8x 5 - Ax*)ti - (14x 4 + 14a; 2 + 6)h + Ax 1 + x 5 + 3a; 3 - 27a;) (5) 



-2x 6 + 2x 4 - 6x 2 )h" (4) 



where h = h(x) is an arbitrary function. The form of the solution is not unique. Because 
h(x) is arbitrary, replacing h by a differential expression in one or more arbitrary functions 
would give a solution too, but a solution with expressions of higher differential order in 
more arbitrary functions without the solution being more general. One naturally seeks a 
general solution which involves only as few as possible arbitrary functions and lowest order 
derivatives of them. But even this requirement does not give a unique solution. One might 
want to minimize the highest order of all parametric functions or, for example, the sum of 
all orders of all parametric functions. 

Furthermore, the occurring function(s) could be scaled to modify the form of the solution, 
for example, to make it denominator free. In the above case replacing h(x) by (x 8 — 2a; 6 + 
7a; 4 — 6x 2 + 9) 3 p(a;) with an arbitrary function p(x) would make the solution polynomial but 
increase its size, i.e. the total number of terms in the coefficients of the parametric functions 
on the right hand sides. 

In the following section the algorithm is given, followed by comments on its character- 
istics. In section |3] and H] we formulate the new algorithm and Euclids algorithm in vector 
notation in order to compare them in detail in section (j5J). We show that both are essentially 



different but also that they compliment each other in the sense that one can give a criterion 
under which circumstances which of both is better suited. Another option is to combine 
both algorithms in a hybrid version. 

Finally, in section ([7]) an application is described which arose from a classification of 
hyperbolic evolutionary vector PDEs. 

2 The Algorithm 

2.1 Outline 

The essence of the algorithm is 

• to partition the ODE into a total derivative and an algebraic remainder, 

• to introduce a new function f r+ i such that the total derivative part is df T+ \/dx and 
thus to write the ODE as a system of 2 equations: one equation defining f r+ i and one 
re- formulating the ODE in terms of f r +i, 

• to use the re-formulated ODE to eliminate and substitute another function and thus 
to arrive again at a single ODE for the same number of functions which in some sense 
is closer to be solvable algebraically than the ODE before. 

These steps are repeated until either the ODE involves only one function and thus is not 
underdetermined anymore, or until one function occurs purely algebraically and thus allow- 
ing the ODE to be solved for that function. The following is a more detailed pseudo code 
description. 

2.2 Pseudocode 

Input • list of functions fx, . . . , f r of x, 
• linear ODE = u in f t (like flTJ) 

Body 

L := {} % L will be a list of substitutions 
s := r 

while (the ODE = u involves at least two fi) and 
(differential orders > V/i) do 

• Factor out 4- from to once as far as possible 
by introducing a new function f s+ i(x) and 
by computing expressions bi. 

s 

= W = / s+1 ' + fi h i + fl 00 (6) 
i=l 

where f s +i,bi are given through: 

s m—i j 

f s+ i = EE#'" H) E(-i)^S (7) 

i=l j=0 k=0 

k = E(-i) (n<+fe) 4r4 } (8) 

k=0 



• if hi = for all i then 

the ODE = u> is exact (apart from aoo), i-e. 
consider new ODE = u := / s+1 + J a 00 <ix 
where is only an abbreviation defined in ([7]) 
else 

o regard / s+ i as a new unknown function, 
o solve (EI) for one function fj of fx, .., f s that 

— has a non- vanishing coefficient bi in ([6]), and 

— is of lowest possible order in oj and thus in ([7]), 

£ = f/^i'+£/A + «oo) (9) 

o use Qj to substitute fj in (J7J) to get a new ODE 

= & ■= —fs+i + 121=1 Sj=o • • • 
for functions f u .., fj_ u f j+1 , .., f s , f s+l 
o update u :— u, s :— s + 1, L := {/j = ..} U L 

end 

if the ODE involves a function / 3 - purely algebraically then 
solve for fj and add it to L : L := {fj = ..} U L 

h(x) := / s (x) % to remember the last introduced function 

P := L % to remember the complete list of substitutions L 

% next substitutions each as stored in the first 

% element of P are performed in the rest of P 

while s > r do 

P := rest(P)\f g= ,„ % as given in first(P) 

s := s — 1 
end 

end % of Body 



Output • I % list of new functions 

• L % the complete list of substitutions 

• P % All initial functions f\, .., f r are either 

% parametric, or are given in P in terms of h 

• = u(x, h(x)) % only if the first while loop 

% terminates due to u not 

% being underdetermined anymore 



2.3 Comments 

Before comparing the algorithm with Euclids algorithm a few comments are necessary. 

• All steps involve only algebra or differentiations with only one exception: if the ho- 
mogeneous part of the ODE is exact then the integral of the inhomogeneous part is 
taken. But this integral does not have to be evaluated, i.e. to be expressed in terms 
of elementary functions. It can stay in a symbolic unevaluated form for the algorithm 
to continue. 

• The coefficients can be arbitrary explicit functions of x, and involve, for example, 
sin or log with the only condition that it must be decidable whether expressions 
involving these functions and their derivatives are zero or not. 



• The transformation of the ODE in each step of the first while loop is reversible, i.e. 
the new ODE is equivalent to the previous one. Therefore the obtained solution is the 
general one. 

• The algorithm terminates. 



We consider how the total sum J2t=i n i of differential orders n, of all functions fa 
changes during execution. After substitution of fj with flU]) in © the new differential 
orders hi of functions fa that occur in the new ODE are: 



The function fj of order rij gets replaced by a new function / s+1 which then occurs 
with same order rij, but the order of all other functions is lowered by at least one 
because we choose j such that rij = min(nfc) (Vfc with bk ^ 0), i.e. we have n, < rij — 1 
also for the fi with 6, ^ 0. Because the ODE has at least two functions, the total sum 
of derivatives is decreasing. The algorithm is therefore finite. 

• The algorithm results in an ODE for a single function iff the differential operators Di 
in (CQ) have a common differential factor. The remaining ODE is of the same order as 
the common factor, i.e. its order can not be higher than the order of the original ODE 



• The algorithm naturally splits into two parts, a part A) establishing a list L of sub- 
stitutions (Q in the first while loop and part B) performing the substitutions in the 
second while loop to obtain an explicit solution. The first part is executed very fast 
(see next section for more details), the second may take longer for higher order ODEs 
because expressions typically grow with each substitution, often exponentially. An 
example is given in the appendix. 

For many applications the list L of substitutions may even be of higher practical value 
than the explicit formulas fi = Fi(x, hi(x), . . . , h r -i(x)), i = 1, . . . , r resulting from 
B). In general, the list L is a much shorter representation of the parametric solution 
of the ODE than the explicit solution itself and thus is more useful as a solution, like 
in the following scenario. 

Let us assume that we have a large algebraic expression in terms of the original 
functions fi,--,f r that has to be evaluated modulo the underdetermined ODE ([I]). 
Instead of replacing the fi directly by their large explicit expressions as given in the 
list P it is usually much better to perform successively the substitutions stored in L 
in the order they were derived, allowing cancellations to happen after each individual 
substitution. Also, if numerical computations are to be done, it is much faster to 
compute the sequence of substitutions than the explicit expressions. 

3 A vector representation 

To perform the iteration process efficiently, i.e. to replace functions by algebraic combina- 
tions of other functions in the ODE (J7J) without needing to perform any differentiations, 
the ODE has to be represented in a form where D := 4- is factored out as far as possible. As 
we will see further down, this is also the appropriate representation for using the Euclidean 



Proof: 




for i = s + 1 
for i < s, bi ^ 
for i < s, bi = 0. 



algorithm to solve the underdetermined ODE. We write the ODE in the form 

r 

= Y. A *f* + a oo(x) (10) 

i 

Ai = DAi + aioix) (11) 
A\ = D^am^x) + . . . + Da i2 (x) + an(x) (12) 

Replacing, for example, f\ = w{x)f\ would only require multiplications a\k '■— wa\k to 
update this representation and no differentiations. 

In this notation a single step in the first while loop of the new method consists of 

• introducing a new function f r+ i through 

/ r+1 = (13) 

i 

giving the ODE the form 

= Df r+1 + ]T a i0 (x)fi + a 00 - (14) 



regarding the defining relation ( I13p for the new function as the new ODE and using 
the old ODE (THj) for a substitution of an 'old' function. For that we choose one a^o of 
the non- vanishing in ffTTl) for which the order rij of the corresponding fj is minimal. 
Thus, performing the substitution 

fj = — — ( Df r+1 + a i0 f t + a 00 ) (15) 

in the definition (TP3|) gives the new ODE 

= ~ Ui + S M - S A > - A > - A ' fe) • (16) 



In vector notation the above two steps change the differential operators (assuming for sim- 
plicity of notation j = 1) according to: 
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(17) 
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V A 
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aoo — > 




-A X 


\a w J 






(18) 



and add (|T5!) to the list of substitutions. Updating components 2 ... r in (|T7|) is done purely 
by multiplications and additions, only updating the first component takes differentiations 
when factoring out D. 

A more conventional algorithm solving underdetermined linear ODEs tries to lower the 
order of the ODE by performing Euclids algorithm until the ODE is algebraic for one 
function and thus can be solved algebraically. In the following section we formulate such an 
Euclidean step to be able to compare it with the new algorithm above. 



4 The Euclid version in vector notation 



For applying the 'right' Euclid algorithm the ODE is given as well in a form with D com- 
pletely factored out: 



= J2 A ifi + a 00<» 

i 

Ai = D ni a ini (x) + . . . + Da a (x) + a i0 (x). 
One iteration step is performed by 



(19) 
(20) 



choosing two functions fa, fj (w.l.o.g. > rij) and introducing a new function f r+ \ 
through 

U = f r+1 - D ni ~ nj (^f) 



(21) 



and performing this substitution in the ODE (|T9|) . 



The substitution is chosen such that the differential order of fa is lowered by at least one. 
This iteration process stops like in the new algorithm when the first function appears purely 
algebraically or when the ODE involves only one function. Therefore, in order to minimize 
the number of steps one would choose fj as one of the lowest order functions and fa as 
one of the lowest order functions from the remaining ones. (From all these choices it is 
recommendable to choose a pair so that multiplication is minimized and Ai, Aj have as few 
as possible terms.) For j — 1, i — 2 the vector notation gives 



( A, \ 

A 2 
A 3 

{ A r J 



— » 



( A r+1 \ 
A 2 
A 3 



( 



At 



A 2 - A\D 
A-, 



ri2—n\ a 2n 2 



(22) 



The update of the 2 nd component in 
as far as possible. 



\ A r ) \ A r J 

aoo -> aoo- (23) 
requires differentiations when D is factored out 



5 Relations between both algorithms 

5.1 Differences 

In this section we want to justify our claim that both algorithms differ significantly from 
each other and that they are not merely variations of one and the same procedure. 

Both algorithms differ conceptually in that Euclids algorithm pairs two differential op- 
erators to get a new differential operator of lower order whereas in the new algorithm the 
relation that defines a new function becomes the ODE and the old ODE is used for substi- 
tuting one of the 'old' functions. 

Another difference is that Euclids algorithm lowers in one step the order of only one 
differential operator whereas the new algorithm lowers the order of all but one differential 
operators. Both types of steps are performed at about the same cost as for the new method 
only the update of the first component in f lT7|) is potentially size increasing and for Euclids 
method it is only the update of the second component of f[2"2"l) . Using Euclids algorithm one 



can of course take one of the lowest order differential operators and decrease the order of 
all others but that would take r — 1 computations, each potentially size increasing. This 
difference does usually not affect the number of necessary iterations (until one function 
occurs purely algebraically), because Euclids method can operate just on two of the lowest 
order operators and ignore all others. So, by applying Euclids method to lower the order 
of the two lowest order operators A, with each other, all other higher order operators keep 
their high order. Moreover, when the algorithm stops because the final equation contains 
one function algebraically, then back substitutions start (part B of the algorithm) which 
even increase the highest orders. 

Differently with the new method where the differential order of all functions but one 
gets lowered in one step. Thus the obtained parametric solution is typically of lower order 
than the solution obtained by Euclids algorithm if this is executed by pairing only lowest 
order derivatives. 

For example, the size of the parametric solution of the ODE 

= x 3 f ax + (x- l)g bx + h 5x , for /, g, h of x (24) 

and the differential order of parametric functions in it depend strongly on the differential 
orders a, b of / and g in the ODE and the method that is used. As a, b increase from 
a = b = ltoa = b = 3 one can see the following trends: 

• In Euclids solution the differential order of parametric functions increases from 5 to 
7 whereas in the solution of the new algorithm the order decreases from 4 to 3 in the 
expression for g, from 4 to 2 in the expression for / and only increases from to 1 in 
the expression for h. 

• The size of expressions in Euclids solution steadily increases and in the solution of the 
new method it decreases. 

For a = b = 1 the solutions have comparable size: 

- When explicit x-dependent factors are not absorbed (for absorbing factors see sec- 
tion 15. 3j) then h is parametric and the rational expressions for /, g have the form (5 
terms)/(2 terms) in Euclids solution and in the solution of the new method / = (4 
terms)/ (2 terms), g — (3 terms) /(2 terms). 

For a = b = 3: 

- When factors are not absorbed then in Euclids solution h is parametric, / = (85 
terms)/ (15 terms), g = 13 terms and in the solution of the new algorithm we have 
/ = 1 term, g = 7 terms, h = 4 terms. 

- When factors are absorbed then Euclids solution gives / = 27 terms, g = 16 terms 
and in the solution of the new algorithms we have / = 1 term, g = 7 terms, h = 4 
terms. 

Not only the solutions of both algorithms may differ significantly but also the number 
of steps to reach them. Starting, for example, with 

= /' + a(x)h {n) = 

to be solved for / = f(x), h = h(x) this takes only one iteration step with the new algorithm 
giving, for example for n = 5, a solution of the form / = (27 terms) /(l term), h — (1 
term)/(l term), whereas it takes n steps for Euclids algorithm giving, for example for n — 5, 
a more spacious solution of the form / = (40 terms) /(l term), h = (2 terms) /(l term). 



If the ODE is inhomogeneous then another difference between both algorithms becomes 
apparent. With non-zero inhomogeneity aoo in the ODE ( TTUj) the transformation ( 1T51) be- 
comes inhomogeneous and changes aoo i n (fl8|) whereas Euclids algorithm does not change 
aoo hi the homogeneous transformation ( )2~Tj) . 

Two ODEs that highlight the duality between both algorithms are 

= /' + f + g'+(ah)W (25) 
= f + f + g' + a(h)W (26) 

where a = a(x) is a given function of x and the ODE has to be solved for /, g, h. As shown 
in table 1 for the ODE fl25|) the new algorithm is to be preferred whereas for ODE ([26]) 
Euclids algorithm gives a shorter solution. 



algorithm 


equation (1251) 


equation (1261) 


new 


f=l term #=22 terms 


f=2 terms g=23 terms 


Euclid 


f=22 terms g=23 terms 


f=2 terms g=3 terms 



Table 1. The size of solutions given by both algorithms (both executed without absorbing 

of factors). 

An explanation of the above behaviour comes from the fact that Euclids method is better 
suited when factors are small, i.e. when the coefficients of the leading derivatives in the 

representation ( ITUl) - (|T2l) effectively cancel each other in a quotient, and less suited when 
these factors are large prime expressions in x. Differently, the new method works best if 
factors ^ are small, i.e. when the coefficients of the algebraic terms effectively cancel each 
other, and less suited when these factors are large prime expressions in x. 

Although both ODEs look similar, they are rather different as (ah)^ has many terms 
if the product rule of differentiation is fully applied and a(h)^ has many terms if is 
factored out as far as possible. Thus both algorithms differ in their suitability for both 
ODEs. 

The above computations and any other tests can easily be performed online (see [7]), 
where access to the computer algebra system REDUCE and to the procedures uode and 
print_uode_solution is provided. 

5.2 Hybridization 

The advantage of having these two, in some sense, complementary algorithms lies in the 
fact that both operate on the same data structure, i.e. both their input and output consists 
of an ODE in the form ( jTOi) - ( 1T21) and both add one substitution to a list of substitutions 
of functions in terms of newly introduced functions. Consequently, both algorithms are 
interchangeable, i.e. one step performed with one algorithm could be followed by another 
step performed with the other algorithm, in order to minimize the size of x-dependent factors 
and thus to minimize growth. A different strategy could be to perform each step with both 
algorithms in parallel, to compare the size of the ODE and/or size of derived substitution 
that both methods give and to choose which one to adopt for this step. This strategy does at 
most double the amount of computation but more likely will lead to savings due to working 
with smaller expressions. 

Any hybrid algorithm performing a mixture of Euclid steps and new algorithm steps is 
still finite because in any such step the sum of differential orders of all functions is decreasing. 

5.3 Absorbing factors 

The following efficiency improving measures work for both methods. A representation ({TO]) 
- ( TL2~1) (and identically (|T9l) - ( )20l) ) of the ODE where D is maximally factored out has the 



advantage that any change of functions fi — > h(x)fi does require only multiplications and is 
done very easily in both algorithms. This freedom of efficiently multiplying functions with 
x-dependent factors can be used for different purposes. 

If in a newly generated ODE all coefficients of a function fi have a non-trivial GCD 
c(x) := GCD(a,i ni , . . . , cijo) 7^ 1 then this can be absorbed into a new function / r+1 and the 
ODE be simplified by performing the substitution 

fi = fr+x/c (27) 

and adding it to the accumulating list of substitutions. Such non-trivial GCDs occur rel- 
atively frequently as the GCD is taken only from the few coefficients of any single one 
operator A4 that changed in an iteration step. 

A different purpose of introducing new functions multiplied with an explicit x-depending 
factor is to avoid a denominator (den). This would arise in the first component of ffTTj) 
and it can be prevented by introducing another new function f r+ 2 and by performing the 
substitution 

fr+i = a±o fr+2- (28) 
Similarly, in the other components of (fTTj) one can scale 

/. = den (*°) f r+l . (29) 

In Euclids algorithm the appearance of denominators in the second component of (122j) 
can be prevented by introducing the scaling 

h = den ( ^) f r+1 . (30) 

\ a lni / 

Another situation where a dominator occurs is at the end of both algorithms when one 
function occurs purely algebraically 

= a i0 fi + A jfj + a oo(x). (31) 

To avoid the denominator a«o one can scale all fj depending on their differential order and 
on their coefficients in Aj. Avoiding this denominator is especially helpful, as fi = ... is 
the last substitution and thus substituted successively in all previous substitutions leading 
easily to a substantial growth of denominators. By avoiding denominators in the above way 
both methods produce denominator free solutions if the ODE is homogeneous. 

5.4 Embedded ODEs 

The case that an underdetermined ODE factorizes, i.e. that the differential operators Ai in 
( flOl) have a non-trivial GCD, or in other words, that the ODE can be written in form of two 
nested ODEs = Q(x, cu(x, fi)) is discovered by both algorithms. The inner ODE = w is 
only determined up to a linear change u> = a(x)u and both algorithms will usually find fi's 
that differ by some a(x). A slight advantage of the new algorithm is that the D° part of 
the ODE is used within the algorithm, thus it is automatically recognized if this vanishes, 
i.e. if the ODE is exact (up to an inhomogeneity) . 



6 ODE-systems 



The introduced methods of solving a single underdetermined ODE (or converting it to an 
ODE for a single function) can be used to convert an ODE system into an equivalent set of 
fully decoupled ODEs, each for a single function. 

After treating any one equation of the original ODE system the original functions in it 
are expressed in terms of fewer functions which are either all free or at most one has to 
satisfy a single ODE. All functions can be replaced in the remaining ODEs. This can go on 
as long as we have equations of at least two functions. When the procedure stops we have 
solved the system or got ODEs, each containing only one function. If more than one ODE 
contain only one and the same single function, then these ODEs form an over-determined 
subsystem which can be treated by a Grobner bases computation and result either in the 
explicit solution for this function or a single ODE for this function of an order not higher 
than the lowest order of the ODEs for this function. This whole procedure terminates with 
either the explicit parametric solution of the original system or a decoupled set of ODEs, 
each ODE for a single function. 

7 An application 

A class of applications where underdetermined linear ODEs occur frequently is the classi- 
fication of hyperbolic evolutionary PDE systems. The aim of such an investigation is to 
find integrable systems of PDEs by determining those systems which have a higher order 
symmetry (see below). More information about the mathematical background is given in 
[6] where a classification of hyperbolic vector PDEs is discussed. 

Let us look at two hyperbolic scalar PDEs for functions u(x, t, r),v(x, t, r) (where x, t are 
the usual independent variables and r is a symmetry parameter). The following ansatz for 
the system and symmetry is generated based on homogeneity considerations. By requiring 
the same homogeneity weights as the potential nonlinear Schrodinger equation we get for 
the system the ansatz 

u tx = a w u 3x + a & u x v 2 x u 2 + a 8 u 2x u x v + a 9 u 2x v x u + aiu 2 x v x + a 3 u 3 x v 2 + a 7 u x v 2x u 

v tx = a 19 v 3x + a 15 u 2 x v x v 2 + a 18 v 2x v x u + a 17 u x v 2x v + a l3 u x v 2 x + a^v^u 2 . (32) 

Assuming the same differential order for the symmetry we get the ansatz 

u T = b w u 3x + v T = b 19 v 3x + (33) 

with right hand sides identical to those of (132]) . only with coefficients hi instead of <Zj. All 
coefficients undetermined functions of the product uvE 

The relations f[3"3"|) are considered to be a symmetry of the system (1321 if the symmetry 
conditions 

d T (u tx ) - d t d x (u T ) = 0, d T (v tx ) - d t d x (v T ) = (34) 

are fulfilled identically in u, v and derivatives of u and v both modulo substitutions based 
on fl32l and fl33l . In performing the differentiations in fl34|) . doing repeatedly substitutions 
(1321 . fl3"3"|) and finally setting all coefficients of different products of powers of derivatives of 
u, v individually to zero gives 27 ODEs with a total of 1334 terms for 26 functions a i; bj of 
z := uv. The length of equations ranges from 2 to 266 terms. In the course of solving this 
system the program Crack performs integrations, substitutions, splittings (i.e. separations 

^or hyperbolic systems the homogeneity weights include negative values, for example here weight(w)=l, 
weight (v)=-l, so that uv has weight zero and therefore we have no limitation on the degree of powers of uv 
and thus all unknown coefficients are arbitrary functions of uv. 



when z occurs only explicitly) and a number of case distinctions. In one of the sub cases the 
resulting conditions can be integrated successively up to the underdetermined linear ODE 

= 36 13 z — 6b 15 z 2 — 2b l7 z 2 + b 17 z — 66152; + 2617. (35) 

Step 1: To start, partial integration gives 

= (3&132: - 6&152: 2 - 2b 17 z 2 + 5b 17 z)' - 36 i3 + 66152 - 3b 17 . (36) 

which can be written as 

= - 36 13 + Qb 15 z - 3b 17 (37) 

by introducing C\(z) through 

ci = 36i3Z — 6bi5Z 2 — 2b 17 z 2 + 5bi 7 z. (38) 

From the functions that occur only algebraically in ( 13?]) (i.e. 613,615,617) the ones that 
have lowest derivatives in ( l35l) are 613,615. Eliminating one of them from (|37j) . say 613 and 
substituting it in ( 138]) gives 

= 2b 17 z 2 - 2b 17 z - c[z + c x (39) 

which is not algebraic in any function yet, but already of first order, so one more step has 
to be performed. 

Step 2: Partial integration of ( 1391) and introduction of 

c 2 = 26i 7 ^ 2 - Cl z (40) 

results in 

= c 2 ' - 6617,2 + 2ci (41) 
which allows to solve for c\ and replace it in (j4"0~]) giving 

= l -c^z - c 2 - b 17 z 2 . (42) 

This condition is purely algebraic for one function, 617, and therefore the algorithm stops. 

Cleanup: We need explicit solutions of (1351) . so what remains to be done are back sub- 
stitutions: (]4"2l provides 

1 / 1 
— c 9 

2z 2 z 

The second substitution expressing 613 in terms of the parametric function c 2 is obtained 
after backward substituting c\ from (14"T!) into (!37|) and 617 from ( 143]) in ( 137|) to get the 
explicit solution consisting of ( 143]) and 

1 // 3 / 2 

fo i3 = gC 2 " - ^c 2 ' + 26152: + — c 2 (44) 

involving the free function c 2 (z). 

Another underdetermined equation resulting in this integrability problem is 

= 36 x 'z - 96 3 V + 36 6 '2: 2 - 46 8 V - 126 8 'z - 126 3 z + 66 6 ^ - 36 8 (45) 

which has the solution 

bx = {Sb^z 2 + 3c^z + 6b 6 z 2 - 5b 8 z - 2c 3 ) /(3z) 
63 = (-46 8 '^ 2 + C32: + 36 6 ,2 2 -3b 8 z-c 3 )/(3z 2 ). 
We finally obtain as a system with higher order symmetries: 

a 1 



.17 = — c 2 - —ci- (43) 



u tx = — ^ (iuv)^v - uv{uv) 2 v x ) 



2uv 3 

a 13 / x 2 
Vtx = —{UV) X V X . 

2uv 
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The following example illustrates the comments made in the last paragraphs of section 12.31 
It shows two possible representations of a solution to an underdetermined ODE. For the 
equation 

(x - l) 3 /i (5) + 3/f } + xf'{ + (1 - x 2 )f[ + fi- (x-2)(x- 3)/£ -xft = 
the generated list L of substitutions is given through 

h = (7289/9/gX 6 - 39400/3/gX 5 + 725945/9/ 8 x 4 - 691667/3/ 8 x 3 + 2665016/9/ 8 x 2 
-463541/3/gX + 40582/g - - 82543/9/ 8 x 5 + 393212/3/ 8 x 4 - 1997777/3/ 8 x 3 
+1404610/ 8 x 2 - 9064364/9/ 8 x + 863648/3/ 8 )/(x 1 2 - 27x x l + 3?>$x l Q - 2551x 9 
+12566x 8 - 43294x 7 + 111667x 6 - 221121x 5 + 332143x 4 - 447477x 3 + 625912x 2 
-422746x + 117948) 

f 7 = (-197/419/^ - 2227/6704/ 6 x 6 + 40655/6704/ 6 x 5 - 282747/6704/ 6 x 4 

+937545/6704/ 6 x 3 - 737203/3352/ 6 x 2 + 236215/1676/ 6 x - 15636/419/ 6 )/ 
(x 5 - 5988/419x 4 + 30423/419x 3 - 64170/419x 2 + 46012/419x - 13152/419) 

/i = (-288/197/7 - 450/197/ 6 )/(x 5 - 2988/197x 4 + 15683/197x 3 - 32571/197x 2 
+18956/197x - 6093/197) 

/ 5 = -9/16/ 6 + 37/32/ix 4 -415/32/!X 3 + 347/8/ix 2 -815/32/ix + 261/32/! 

h = -4/9/^ - 17/9/iX 3 + 124/9/iX 2 - 176/9/ x x + 11/j 

h = -1/4/1 + 5/4/ix 2 - 17/4/ix - ll/4/i 

h = fs + 2/ix + /! 

which is a much shorter representation of the solution than the 10 page explicit form which 
results from substituting / 6 , / 7 , fx, f 5 , / 4 , / 3 in this order into each other and takes the form 
fx =(42 terms)/(25 terms), / 2 =(304 terms)/(61 terms) involving up to 29-digit integers. 
The list of 7 substitutions is not only shorter than the list of 2 substitutions for / 1; / 2 , it 
also is much faster to derive and more useful if a differential expression is to be simplified 
modulo the solution of the above ODE by substituting / 2 , fs, /s, fx, fi, fe in this order. 

The difference in size of both solution representations can be arbitrarily amplified by 
having an input ODE of higher order and higher degree polynomials as coefficients. 

9 Summary 

We present an algorithm for the solution of underdetermined linear ODE that is compatible 
but structurally different from the (right) Euclidean algorithm. Because both algorithms 
operate on the same data structure and because both complement each other in the sense 
that each one is most efficient for ODEs of different form, the combination of both algorithms 
is superior to each individual one. 
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